clc;
clear;

animate = 0;

n = [3 4 5 6 7 8 12 16 32 64 128];
exact = @(x) sqrt(x.*x + 1);

figure;
hold on;
x = 0:0.01:1;
plot(x, exact(x), '-r');

videoObj = VideoWriter('Euler_animate.mp4', 'MPEG-4');
videoObj.FrameRate = 2;

if animate
    open(videoObj);
end

for i = 1:length(n)
    h = 1/n(i);
    x = 0:h:1;
    y = [];
    y(1) = 1;
    for j = 1:length(x) - 1
        y(j + 1) = y(j) + h * (x(j) / y(j));
    end
    animate_plot(x, y, n, animate, videoObj);
end

if animate
    close(videoObj);
    close(gcf);
end

function animate_plot(x, y, h, animate, videoObj)

if ~animate
    hold on;
    xmin = min(x);
    xmax = max(x);
    ymin = min(y);
    ymax = max(y);
    pad = 0.05;
    xmin = xmin - pad*(xmax-xmin);
    xmax = xmax + pad*(xmax-xmin);
    ymin = ymin - pad*(ymax-ymin);
    ymax = ymax + pad*(ymax-ymin);
    axis([xmin xmax ymin ymax]);
    grid on;
    title(['step: ', num2str(h)]);
    xlabel('x');
    ylabel('y(x)');
    plot(x, y, '-b');
    legend({'Exact', ['Euler:', num2str(h)]})
else
    xmin = min(x);
    xmax = max(x);
    ymin = min(y);
    ymax = max(y);
    pad = 0.05;
    xmin = xmin - pad*(xmax-xmin);
    xmax = xmax + pad*(xmax-xmin);
    ymin = ymin - pad*(ymax-ymin);
    ymax = ymax + pad*(ymax-ymin);
    axis([xmin xmax ymin ymax]);
    grid on;
    title(['step: ', num2str(h)]);
    xlabel('x');
    ylabel('y(x)');
    plot(x, y, '-b');
    legend({'Exact', ['Euler:', num2str(h)]});
    frame = getframe(gcf);
    writeVideo(videoObj, frame);
end

end